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A  model  of  homogeneneous  nucleation  based  on  first  pronciples  of  the  kinetic  theory 
of  gases,  and  developed  for  the  direct  simulation  Monte  Carlo  method,  is  presented.  Key 
parameters  of  the  model,  such  as  the  binding  energy  and  the  heat  capacity  as  function  of 
cluster  size,  are  given  for  argon  and  water.  For  water,  two  datasets  are  used,  one  based  on 
quantum  calculations  of  heat  capacity  and  binding  energy  available  in  the  literature,  and 
the  other  based  on  the  present  molecular  dynamics  computations.  The  model  is  analyzed 
using  the  thermal  bath  relaxation  problem,  and  applied  to  compute  the  water  dimer  teminal 
mole  fractions  in  circular  orifice  expansion. 

I.  Introduction 

Two  different  approaches  for  modeling  the  condensation  in  rapidly  expanding  plumes  have  been  reported 
in  the  literature.  The  first  approach,  known  as  the  classical  approach,  takes  its  starting  point  from  the 
classical  nucleation  theory  (CNT)  which  is  based  on  equilibrium  thermodynamics.1,2  Alternatively,  one  can 
treat  nucleation  as  the  process  of  kinetic  cluster  aggregation.3  A  discussion  of  the  assumptions  and  difficulties 
with  using  CNT  for  predicting  condensation  in  supersonic  expansions  can  be  found  in  earlier  papers. 4-11 

A  promising  direction  in  modeling  the  coupled  condensation  flow  is  the  use  of  a  kinetic  particle  simulation 
method,  direct  simulation  Monte  Carlo  (DSMC),12  which  is  applicable  in  a  wide  range  of  flow  regimes 
from  free  molecular  to  near  continuum.  The  DSMC  method  has  been  used  to  study  the  process  of  cluster 
formation  and  evolution  for  a  number  of  years. 13-15  However,  the  gas  flow  in  the  earlier  studies  was  uniform, 
the  considered  cluster  size  range  was  very  narrow  (up  to  25  monomers  in  a  cluster)  and  the  examined 
reaction  types  were  unrealistically  limited  to  elastic  collisions,  cluster  and  monomer  sticking  to  clusters, 
and  evaporation  of  monomers  from  clusters.  More  recently,  the  DSMC  method  has  been  extensively  and 
successfully  applied  to  modeling  the  processes  of  cluster  formation  and  evolution  in  supersonic  jets  by  Levin 
et  al  (see,  for  example,  Refs.  9,  16, 17).  The  model  initially  was  based  on  CNT,  with  the  new  clusters 
being  formed  at  the  critical  size.  Further  work  of  these  authors18  extended  the  kinetic  dimer  formation 
approach  of  Ref.  19  (which  assumed  that  a  ternary  collision  always  results  in  a  dimer  formation),  to  include 
molecular  dynamic  (MD)  simulations  for  obtaining  information  on  the  probability  of  dimer  formation  in 
ternary  collisions.  The  work20  used  a  temperature-dependent  probability  of  formation  of  argon  dimers.  In 
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addition,  these  authors  used  MD  to  investigate  the  kinetics  of  H20  dimer  formation  and  showed  evidence  for 
a  bimolecular  mechanism21  and  key  predictions  were  made  (using  MD)  on  many  needed  kinetic  parameters 
for  argon16  and  H2022  cluster  formation.  None  of  the  previous  papers  as  yet  have  combined  the  kinetic 
nucleation  mechanism  and  the  kinetic  evaporation  mechanism  together. 

A  previous  paper23  introduced  the  present  method,  where  the  DSMC  approach  for  modeling  of  homo¬ 
geneous  nucleation  in  rapidly  expanding  plumes  is  extended  to  include  a  number  of  new  features.  Most 
importantly,  a  kinetic  RRK  model24  is  implemented  to  characterize  the  cluster  evaporation  rates.  Then, 
an  energy  dependent  collision  procedure  similar  to  the  recombination  reaction  model  of  Ref.  25  is  used  for 
the  collision  complex  formation.  An  empirical  parameter  is  used  for  the  inelastic  collision  number  in  the 
cluster-monomer  collisions.  For  dimers,  this  parameter  is  calibrated  through  the  comparison  of  the  computed 
nucleation  rates  and  equilibrium  constants  in  thermal  bath  with  available  theoretical  and  experimental  data 
for  argon  and  water.  The  present  paper  briefly  recapitulates  the  method  used  and  presents  a  number  of  up¬ 
dates  to  the  argon  and  water  parameters  used.  The  importance  of  physically  realistic  values  for  parameters 
is  emphasized  and  new  MD  results  for  sample  parameter  values  are  presented,  along  with  a  discussion  of  the 
sensitivity  of  results  to  these  values.  Finally,  a  new  validation  case  is  presented  by  comparison  to  published 
data  for  H20  dimer  measurements. 


II.  Numerical  Model 

The  condensation  model  used  in  this  work  is  formulated  for  the  DSMC  method  and  is  based  on  first 
principles  of  the  kinetic  theory.  The  considered  processes  of  cluster  nucleation  and  evolution  are  modeled  at 
the  microscopic  level.  First  principles  theory  are  used  to  define  the  main  processes  of  homogeneous  conden¬ 
sation,  where  all  collision,  nucleation,  and  evaporation  events  depend  on  instantaneous  energies  of  colliding 
partners,  and  not  cell  temperature  or  other  macroscopic  quantities.  The  processes  that  are  included  in  the 
model  and  outlined  below  are  (i)  formation  of  collision  complexes  through  the  binary  collisions  of  cluster¬ 
forming  monomer  species  and  creation  of  dimers  through  the  collision  stabilization  of  collision  complexes,  (ii) 
elastic  monomer-cluster  collisions  that  change  the  translational  and  internal  energies  of  colliding  particles, 
(iii)  inelastic  monomer-cluster  collisions  that  result  in  monomer  sticking,  (iv)  cluster-cluster  coalescence,  (v) 
evaporation  of  monomers  from  clusters.  The  description  of  each  of  these  processes  is  given  below,  with  all 
necessary  constants  and  variables  specified  for  two  gases  considered  in  this  work,  argon  and  water;  more 
detail  on  the  algorithm  may  be  found  in  Ref.  23. 

Dimer  Formation  One  of  the  important  assumptions  of  the  present  model  is  that  all  pairs  of  colliding 
particles  create  collision  complexes.  A  collision  complex  is  a  pair  of  monomers  that  have  collided,  and  may 
have  the  conditions  necessary  to  form  a  dimer  if  struck  by  a  third  particle  during  its  lifetime.  The  collision 
complex  lifetime,  ti ,  is  assumed  to  be  dependent  on  the  type  of  monomers  and  their  relative  collision  velocity, 
with  the  functional  dependence  given  by  the  well  known  Bunker’s  expression26 

ti  =  1.5ct0/li5  (1) 

where  op  and  eo  are  the  potential  depth  and  separation  distance  parameters  of  the  Lennard-Jones  (LJ) 
potential,  [i  is  the  reduced  mass  of  the  colliding  particles,  and  E  is  their  relative  translational  energy.  The 
values  of  op  and  eo  used  in  this  work  are  3.2  x  10~10  m  and  7.94  x  10-21  J  for  water,  and  3.405  x  10“10  m  and 
1.654  x  10-21  J  for  argon.  The  argon  LJ  parameters  are  well  known,  but  H20  values  are  roughly  estimated 
by  reproducing  the  known  viscosity  of  water  vapor  at  273  K  and  for  a  small  range  of  temperatures  using  a 
LJ  potential  model. 

The  process  of  interaction  of  collision  complexes  with  surrounding  gas  particles  is  modeled  using  the 
majorant  frequency  scheme27  with  the  assumption  that  the  collision  complex  -  third  particle  interactions 
are  governed  by  the  Variable  Hard  Sphere  (VHS)  interaction  model.28  The  VHS  viscosity-temperature 
exponent  to  of  a  collision  complex  was  assumed  to  be  that  of  the  comprising  monomers.  For  argon  atoms, 
the  VHS  parameters  taken  from  Ref.  12  were  assumed,  dref  =  4.17  A  and  u>  =  0.81.  For  water  molecules, 
dref  =  6.2  A  and  ui  =  1  were  assumed,  based  on  reproducing  the  known  viscosity  of  water  vapor  at  273  K 
and  for  a  small  range  of  temperatures  using  the  VHS  potential  model.  Note  that  the  VHS  collision  diameters 
are  used  to  compute  collision  frequency  in  the  cell,  while  the  LJ  diameters  are  only  used  to  compute  the 
diameter  and  lifetime  of  a  collision-pair  complex  for  purposes  of  determining  the  probability  of  a  three-body 
collision. 
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For  consistency  with  the  collision  complex  lifetime  determination,  an  expression  for  the  diameter  d  of  the 
collision  complex  recommended  in  Ref.  26  was  used, 

d  =  3^a0{e0/E)i.  (2) 

If  there  is  a  physical  collision  between  a  collision  complex  and  a  third  particle,  there  is  a  possibility  of 
forming  a  stable  dimer  as  a  result  of  such  a  collision.  Generally,  the  probability  of  the  formation  of  a  stable 
dimer  (dimer  stabilization)  depends  on  the  colliding  species  and  the  energies  -  both  translational  and  internal 
-  of  the  colliding  particles.  In  this  work,  constant  stabilization  probabilities  of  0.25  for  Ar18  and  0.7  for  H2O 
were  assumed,  which  seem  reasonable  for  the  range  of  temperatures  under  consideration. 

Dimer  creation  through  the  collisional  stabilization  of  collision  complexes  is  modeled  as  a  two-step  process, 
L  +  M  — >  {LAI),  {LAI)  +  K  — ►  LAI  +  K.  Here,  L  and  M  are  monomers,  {LAI)  is  the  collision  complex,  and 
K  is  the  third  particle. 

Reflective  collisions  of  monomers  and  clusters  The  collisions  between  monomers  and  clusters  is 
one  of  the  key  processes  that  determine  the  nucleation  rate.  The  reason  for  this  is  strong  dependence  of  the 
evaporation  rate  on  the  cluster  internal  energy.  Since  the  monomers  are  dominant  in  the  flows  considered 
here,  the  cluster  internal  energy  is  mostly  governed  by  its  relaxation  through  cluster-monomer  collisions. 
In  this  work,  a  hard  sphere  model  is  assumed  for  cluster-monomer  collisions,  with  the  cluster  diameter 
determined  from  Eqn.  2  for  dimers,  and  for  larger  clusters  from  an  empirical  correlation  used  extensively  in 
the  past  (see,  for  example,  Ref.  8), 

d=  2-  {A-  *s  +  £),  (3) 

where  A  and  B  are  species-dependent  constants,  and  i  is  the  number  of  monomers  in  the  cluster.  In  this 
work,  the  values  of  A  and  B  were  2.3  x  10-10  m  and  3.4  x  10~10  m  for  argon,16  and  1.9  x  10~10  m  and 
2.4  x  10_1°  m  for  water.29 

For  the  energy  transfer  between  the  relative  translational  and  internal  modes  of  the  cluster  and  monomer, 
the  Larsen-Borgnakke  model30  is  used,  and  a  parameter  Z  is  introduced  that  has  a  meaning  of  the  internal 
energy  relaxation  number.  The  energy  transfer  between  all  energy  modes  of  the  cluster-monomer  pair  occurs 
with  a  probability  Z _1,  and  an  elastic  collision  with  no  internal  energy  exchange  occurs  with  the  additional 
probability  1  —  Z^1 .  For  argon,  temperature  dependent  values  of  Z{T)  were  used,  obtained  through  the  linear 
interpolation  between  the  values  given  in  Table  1.  Similar  to  temperature  dependent  collision  numbers  for 
rotational  and  vibrational  relaxation  of  molecules  widely  used  in  the  DSMC  method,  T  was  the  cell-based 
translational  temperature.  Such  a  temperature  dependence  allows  good  agreement  of  the  DSMC  rates  for 
dimer  nucleation  and  dissociation  with  rates  available  in  the  literature.  For  water  condensation,  a  constant 
value  of  Z  =  10  was  used,  which  corresponds  to  that  for  rotational  and  vibrational  relaxation  of  water 
molecules  at  romm  temperature  obtained  in  molecular  dynamics  studies.31,32 


T,  K 

0.0 

100.0 

200.0 

300.0 

400.0 

500.0 

Z"1 

0.25 

0.13 

0.08 

0.06 

0.046 

0.04 

Table  1.  Inelastic  collision  number  as  function  of  temperature  for  argon. 


Sticking  of  MONOMERS  and  clusters  When  a  monomer  collides  with  a  cluster,  sticking  of  the  monomer 
to  a  cluster  surface  is  possible,  in  addition  to  a  reflective  collision  described  in  the  previous  section.  For  small 
clusters,  monomer  sticking  is  the  main  process  that  governs  the  evolution  of  the  droplet  size  distribution.33 
For  water  molecules,  an  empirical  dependence  of  the  sticking  probability  on  the  species  radius  and  mass, 
given  in  Ref.  34,  is  used.  This  dependence  reduces  to 

_  d2n  (  mn 
{dn  +  d\)2  \mn  +  mi 


where  indices  n  and  1  refer  to  the  cluster  and  monomer,  respectively. 

For  argon,  the  sticking  probability  of  monomers  on  clusters  given  in  Table  2  is  used 
Similar  to  monomer-cluster  collisions,  the  outcome  for  cluster-cluster  collisions  is  assumed  to  be  either 
coalescence  or  elastic  interaction.  The  probability  of  sticking  was  assumed  to  be  unity  both  for  argon  and 
water  in  cluster-cluster  collisions. 
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Cluster  size 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

Probability 

0.06 

0.075 

0.1 

0.16 

0.3 

0.5 

0.67 

0.75 

0.8 

0.833 

0.857 

0.876 

0.891 

0.9 

Table  2.  Sticking  probability  as  function  of  cluster  size  for  argon. 


Cluster  evaporation  Following  Ref.  35,  RRK  theory  is  used  to  model  the  evaporation  process.  The 
evaporation  rate  ke  is  calculated  as 

*,  =  ,i (5) 
\  i^int  J 

Here,  n  is  the  number  of  monomers  in  the  cluster,  v  is  the  vibration  frequency,  Ns  is  the  number  of  surface 
atoms,  and  Eint  is  the  cluster  internal  energy.  For  dimers,  the  exponent  3n  —  7  is  replaced  with  1.  The 
number  of  surface  atoms  is  n  for  N  <  5,  n  —  1  for  4  <  n  <  7,  and  (367t )1/3(n1//3  —  l)2  for  n  >  6.  The 
vibration  frequency  was  taken  to  be  2.68  x  1012  s-1  for  water  clusters,36  and  1012  s_1  for  argon  clusters.36 

It  is  important  to  use  reasonable  values  for  the  evaporation  energy  of  a  monomer  off  a  cluster,  as  well  as 
the  number  of  cluster  internal  degrees  of  freedom,  as  a  function  of  cluster  size.  In  this  work,  two  data  sets 
were  used  for  water  heat  capacities  and  binding  energies.  The  first  one  is  described  in  this  section;  hereafter 
it  is  referred  to  as  Dataset  1.  The  second  data  set  is  based  on  water  dimer  and  trirner  heat  capacities 
calculated  with  the  approach  presented  in  the  next  section. 

The  evaporation  energy  of  a  monomer  from  an  n-cluster  is  given  by  Eb(n)  —  Eb(n  —  1),  where  Eb  is  the 
binding  (evaporation)  energy  of  the  n-cluster.  The  Eb  values  for  water  clusters  (Ref.  37)  were  taken  from 
Ref.  38  for  n  =  2  —  8  (and  subtracting  the  zero-point  energies,  taken  from  ref  Ref.  39),  and  using  smoothed 
values  from  Refs.  40,41  for  n  =  9  —  13.  These  baseline  values  of  the  evaporation  energy  for  water  clusters 
are  used  in  (Dataset  1).  Evaporation  energies  were  taken  from  Ref.  42  for  argon  clusters.  The  number  of 
cluster  internal  degrees  of  freedom  is  calculated  from  the  expression  for  the  average  internal  energy  (E)  of  a 
cluster  of  a  size  n 

tint  o 

(E)  =  ^-kT  =  nCvT--kT 
as 

ent  =  n2^  -  3 
k 

where  Cv  is  the  cluster  heat  capacity.  For  water,  the  values  of  Cv  where  adapted  from  Refs.  43,44  (Dataset 
1),  while  for  argon,  they  where  assumed  to  approximate  an  expression  £ mt  =  2(3 n  —  r))  +  e,  where  =  5 
and  e  =  2  for  n  =  2  and  77  =  6  and  e  =  3  otherwise.45  Note  that  for  water  dimers,  the  number  of  degrees  of 
freedom  was  assumed  t  be  a  function  of  internal  energy  (or  effective  temperature),  decreasing  linearly  from 
its  listed  value  at  200  K  to  3  at  0  K.  Both  the  heat  capacity  and  evaporation  energies  are  listed  in  Table  3. 
For  the  cluster  sizes  larger  than  given  below,  the  values  for  the  maximum  listed  sizes  are  used. 

The  first-principles  condensation  model  described  here  was  implemented  in  the  DSMC  code  SMILE.46 

III.  Calculation  of  Water  Dimer  and  Trimer  Heat  Capacities 

Water  clusters  have  been  subject  to  several  experimental  and  theoretical  studies.  Tsai  and  Jordan47,48 
studied  the  phase  transition  of  water  octomers.  The  single  histogram  method49  and  jump-walking  proce¬ 
dure50  were  used  in  those  studies.  The  former  allows  one  to  perform  a  single  Monte  Carlo  simulation  at 
some  specified  temperature  and  extrapolate  the  results  to  other  nearby  temperatures.51  The  latter  carries 
out  calculations  for  a  range  of  temperature,  starting  at  a  ‘high1  temperature  and  stepping  to  lower  temper¬ 
atures.  The  initial  calculation  at  temperature  T)  is  carried  out  using  standard  Metropolis  sampling.  Then, 
at  each  successive  temperature  T),  sampling  is  done  by  occasionally  jumping  to  a  distribution  generated  at 
the  preceding  temperature  T)_  1.  TIP3P52  and  TIP4P53  potentials  were  used  in  those  Monte  Carlo  simula¬ 
tions.  Plots  of  the  potential  energy  and  heat  capacity  as  a  function  of  the  temperature  were  generated.  A 
Molecular  Dynamics  study  of  (#20)n=2,3,4,6,8  clusters  was  conducted  by  Guvenc  and  Anderson.54  Melting 
temperature  was  plotted  as  a  function  of  cluster  size,  and  the  ones  for  the  dimer  and  the  tetramer  closely 
ressemble  the  bulk  melting  temperature,  while  those  for  the  other  sizes  were  considerably  lower.  Pedulla  and 
Jordan55  studied  the  melting  behavior  of  (#20)6  and  (H20)$  water  clusters.  The  location  and  sharpness 
of  the  melting  transition  were  investigated  for  different  potential  models.  It  was  found  that  the  position 
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Cluster  size 

H20  Eb,  J 

h2oc„,  j/k 

Ar  Eb,  J 

Ar  Cv,  J/K 

2 

2.455E-20 

3.726E-23 

1.98E-21 

2.41E-23 

3 

5.352E-20 

4.830E-23 

3.96E-21 

2.76E-23 

4 

6.325E-20 

5.748E-23 

5.94E-21 

3.10E-23 

5 

4.726E-20 

6.555E-23 

6.08E-21 

3.31E-23 

6 

4.587E-20 

7.245E-23 

6.89E-21 

3.44E-23 

7 

5.560E-20 

7.817E-23 

7.49E-21 

3.54E-23 

8 

6.950E-20 

8.277E-23 

6.36E-21 

3.62E-23 

9 

5.560E-20 

8.648E-23 

8.29E-21 

3.68E-23 

10 

5.699E-20 

8.970E-23 

8.25E-21 

3.72E-23 

11 

5.838E-20 

9.227E-23 

8.27E-21 

3.76E-23 

12 

5.977E-20 

9.467E-23 

9.86E-21 

3.79E-23 

13 

6.116E-20 

9.654E-23 

12.2E-21 

3.82E-23 

Table  3.  The  water  and  argon  cluster  heat  capacities  and  evaporation  energies  per  monomer. 


of  the  peak  in  the  heat  capacity  curve  was  quite  sensitive  to  the  specific  model  potential.  The  results  ob¬ 
tained  when  using  canonical  and  microcanical  ensembles  for  Monte  Carlo  simulations  applied  to  study  water 
tetramer  and  octamer  were  compared.56  A  TIP5P  potential  model57  was  implemented.  The  microcanonical 
heat  capacity  curve  was  found  to  have  a  sharper  peak  than  the  canonical  one,  but  the  locations  of  the  peaks 
were  found  to  be  the  same.  The  phase  change  of  water  octamer  clusters  was  investigated58  by  performing 
Molecular  Dynamics  simulations  with  a  SPC/F2  potential  model.59  The  phase  transition  is  found  to  occur 
around  125  K. 


A.  Simple  Point  Charge  Water  Model 

The  interactions  among  water  molecules  are  dominated  by  dipole  interactions.  One  effective  way  to  describe 
such  interactions  is  to  consider  three  point  charges,  one  on  each  atom.  In  the  SPC  model,60  the  water 
molecule  is  modeled  to  have  three  centers  of  concentrated  charge:  a  positive  charge  on  each  of  the  two  H 
atoms  and  a  negative  charge  on  the  O  atom.  The  assumption  that  there  are  point  charges  is  an  approximation 
that  leads  to  an  incorrect  value  for  the  permanent  dipole  moment  of  the  water  molecule.  To  correct  this, 
the  H-O-H  bond  angle  is  changed  from  the  true  value  of  104.45  to  109.47  degrees  in  the  SPC  model.  As 
a  consequence  of  the  charge  concentration  and  the  widened  V-shaped  bond  angle,  the  permanent  dipole 
moment  of  the  SPC-model  water  molecule  has  a  value  close  to  that  measured  in  experiment  of  1.85  D. 

In  summary,  the  SPC  model  consists  of  a  triangular  water  model  with  an  OH  distance  of  1  A  (compared 
to  the  true  bond  length  of  0.9584  A),  with  point  charges  on  the  oxygen  and  hydrogen  positions  of  -0.82  and 
+0.41  e  (electronic  charge  units),  respectively.  The  corresponding  potential  is  a  combination  of  Lennard- 
Jones  interactions  between  the  oxygen  atoms  of  each  water  molecule, 


C/LJ(r)  =4eLJ 


(6) 


with  the  parameters  of  ctlj  =  3.166  A,  eLJ 


0.65  kJ/mol  and  the  Coulomb  potential  between  all  the  atoms, 


£|C°ul°mb(r_) 


47re07 


(7) 


where  qi ,  qj  are  the  charges  of  O  or  H  atoms  and  eo  is  the  permittivity  of  free  space.22 

B.  Monte  Carlo  Canonical  Ensemble  Simulation 

For  cluster-monomer  interactions,  the  initial  separations  between  molecules  inside  a  cluster  and  their  ve¬ 
locities  have  a  great  effect  on  the  calculation  of  the  potential  energy.  Therefore  a  high  number  of  initial 
configurations  (between  400  and  4000)  is  chosen  to  reduce  the  statistical  error.  In  order  to  prepare  the 
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initial  position  configurations  for  clusters,  the  Monte  Carlo  canonical  ensemble  simulation  method  is  used.61 
In  this  approach,  a  sphere  domain  is  setup  with  a  known  radius,  Ri,  which  depends  on  the  cluster  size,  i. 
The  initial  molecules  (for  example,  2  molecules  for  a  dimer)  then  are  randomly  put  in  the  spherical  domain 
and  one  of  the  molecules  is  moved  a  small  distance  from  its  original  point  as  shown  in  Figure  1.  The  system 
potential  energy  is  calculated  based  on  the  SPC  model  and  designated  as,  UQ.  The  system  is  allowed  to 
evolve  for  a  large  number  of  steps  and  for  each  time  step,  one  of  the  molecules  is  randomly  moved  in  the 
sphere  and  the  new  energy  Un  is  calculated.  The  move  is  accepted  with  probability, 

Pa  =  min (1,  exp (-[[/„  -  U0))/kT)  (8) 

where  T  is  the  temperature  and  k  is  Boltzmann’s  constant.  After  10,000  steps,  the  system  is  assumed  to  be 
in  equilibrium  and  an  accepted  configuration  is  recorded  every  100  steps.  The  so-called  ’baby  steps’  ensure 
a  higher  number  of  configurations  are  accepted.  The  molecule  is  also  rotated  with  a  small  angle  change 
compared  to  the  angle  at  the  previous  step  as  shown  on  Figure  2.  The  system  is  run  for  1,000,000  timesteps. 
Six  independents  Monte  Carlo  cycles  were  computed  at  every  temperature. 


Figure  1.  Water  trimer  plot  showing  the  outside  boundary  sphere  (radius  =  2.8  A)  and  an  imaginary  cube 
(length  =  0.28  A)  inside  which  ’baby  translations’  are  allowed. 


Figure  2.  The  water  molecule  is  rotated  through  a  small  angle  equal  to  0.27T. 


C.  Potential  Energy  and  Heat  Capacity  Calculation 

The  potential  energy  of  the  system  is  calculated  by  summing  the  SPC  potential  due  to  each  molecule  in  the 
cluster.  Then  it  is  normalized  by  dividing  by  the  number  of  molecules  in  the  cluster.  The  constant  volume 
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heat  capacity  is  computed  as 


Cv 


<  u2  >  -  <  u  >2 

3  RT2 


+  37? 


(9) 


where  R  is  the  ideal  gas  constant.  The  standard  deviations  are  based  on  the  6  independent  cycles  at  each 
temperature. 


D.  Results  of  Molecular  Dynamics  Calculations 

Figure  3  shows  the  potential  energy  of  a  water  molecule  in  a  dimer  as  a  function  of  the  ensemble  temperature. 
It  can  be  seen  that  the  potential  energy  increases  with  the  temperaure,  almost  with  a  linear  trend.  Figure  4 
shows  the  constant  volume  heat  capacity  as  a  function  of  the  temperature.  Even  though  the  value  of  Cv 
fluctuates  a  lot  with  temperature,  the  general  trend  that  can  be  seen  from  this  plot  is  that  there  is  a  peak 
around  220  K.  This  point  where  the  heat  capacity  reaches  reaches  its  maximum  is  known  as  the  melting 
point.  This  value  is  relatively  close  to  the  one  found  by  Guvenc  and  Anderson,54  275  K.  The  heat  capacity 
fluctuates  around  10  cal/mol/K,  while  the  one  for  liquid  water  is  about  18  cal/mol/K  at  300  K. 

The  results  for  the  water  trimer  are  shown  in  Figures  5  and  6.  The  same  linear  trend  as  for  the  dimer 
can  be  seen  again  in  the  case  of  the  trimer.  The  values  of  the  potential  energy  are  about  double  those  of 
the  dimer.  The  trend  for  the  heat  capacity  of  the  trimer  is  also  similar  to  the  one  of  the  dimer.  The  heat 
capacity  fluctuates  around  10  cal/mol/K  and  seem  to  have  a  melting  point  around  240  K. 

In  this  work,  the  dimer  and  trimer  heat  capacities  of  10.2  cal/mol/K  and  12  cal/mol/K  were  used  in 
Dataset  2.  The  binding  energy  of  3e  —  23  J  was  used  for  dimer  in  order  to  reproduce  theoretical  equilibrium 
constants  (see  below);  for  larger  clusters,  the  binding  values  from  Table  3  were  taken. 


Figure  3.  Potential  energy  of  a  water  molecule  inside  a  dimer  for  different  cluster  temperatures.  The  standard 
deviations  are  computed  from  six  2-million  cycles  of  calculations. 


IV.  Thermal  bath  relaxation 

Inelastic  cross  sections  for  monomer-monomer  and  monomer-cluster  collision  processes  are  necessary  for 
detailed  validation  of  a  kinetic  condensation  model.  These  cross  sections,  generally  functions  of  the  energy 
states,  both  translational  and  internal,  of  pre-  and  post-collisional  particles,  are  not  available  for  most  gas  and 
temperature  conditions  of  interest.  Contrary  to  the  energy  dependent  cross  section,  the  integral  temperature 
dependent  rates  for  such  collisions  at  conditions  close  to  equilibrium  are  available  in  the  literature.  Therefore, 
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Figure  4.  Constant  volume  heat  capacity  of  a  water  molecule  inside  a  dimer  for  different  cluster  temperatures 
(heat  capacity  per  molecule). 


Figure  5.  Potential  energy  of  a  water  molecule  inside  a  trimer  for  different  cluster  temperatures.  The  standard 
deviations  are  computed  from  six  1-million  cycles  calculations. 


one  of  the  key  indicators  of  the  accuracy  and  reliability  of  a  condensation  model  is  its  ability  to  produce 
realistic  rates  of  evaporation  and  nucleation  at  equilibrium.  Although  matching  the  rates  generally  does  not 
guarantee  correct  behavior  in  nonequilibrium,  it  still  is  a  necessary  condition  for  a  model  to  satisfy. 

In  this  work,  thermal  bath  relaxation  of  pure  argon  and  pure  water  are  examined  at  different  temperature 
conditions,  and  the  dimer  formation  rates  for  argon  and  equilibrium  constants  for  the  formation  of  dimers 
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Figure  6.  Constant  volume  heat  capacity  of  a  water  molecule  inside  a  trimer  for  different  cluster  temperatures 
(heat  capacity  per  molecule). 


in  argon  and  water  are  calculated  and  compared  to  the  published  results.43, 62-65  In  all  thermal  bath  results, 
one  million  simulated  particles  were  used,  and  the  run  proceeded  until  the  steady  state  is  reached,  after 
which  the  results  were  sampled  for  20  thousand  timesteps.  The  number  density  was  5  x  1023  molec/m3  for 
argon  and  2  x  1023  molec/m3  for  water.  The  timestep  of  2.5  x  10-11  s  was  selected  so  that  the  number  of 
collisions  per  molecule  is  much  smaller  than  unity,  and  the  results  are  independent  on  the  timestep. 

The  computed  dimer  formation  rates  krec  for  argon  are  presented  in  Fig.  7  and  compared  with  the  stable 
dimer  formation  rate  of  Ref.  62,  where  they  were  calculated  using  classical  trajectories,  and  the  following 
expression  was  proposed, 

krec  =  10.15T-0,278  exp  {— 0.0031T}  . 

Generally,  the  present  dimer  formation  rates  are  in  good  agreement  with  the  classical  trajectory  calculations, 
with  the  maximum  difference  approaching  20%  for  higher  temperatures.  Note  that  the  computed  rate  has 
somewhat  different  slope  than  that  of  Eqn.  IV.  A  number  of  factors  could  be  affecting  the  slope,  among 
them  are  the  energy  dependence  of  the  stabilization  probability  and  dimer  heat  capacity,  which  were  not 
included  in  the  present  model. 

The  computed  equilibrium  constant,  Keq,  that  is  the  ratio  of  the  dimer  dissociation  to  the  dimer  formation 
rate,  is  given  in  Fig.  8.  It  is  compared  to  the  theoretical  results  of  Ref.  63,  where  a  number  of  approximate 
classical  and  quantum  methods  are  compared  with  exact  numerical  calculations,  and  also  experimental 
results  of  Ref.  64  and  theoretical  results  of  Ref.  65.  Note  that  while  the  results  for  different  models  and 
interaction  potentials  where  found  to  be  widely  different  in  Ref.  63,  there  was  a  good  agreement  between 
analogous  quantum  and  classical  calculations.  Figure  8  shows  that  for  the  entire  temperature  range  there  is 
an  excellent  agreement  between  the  present  model  and  the  theoretical  and  experimental  values.  The  reasona 
for  such  a  good  agreement  is  the  appropriate  selection  of  the  temperature  dependence  of  the  inalastic  collision 
number  Z. 

This  parameter,  which  s  in  effect  the  inverse  probability  of  the  energy  transfer  between  the  internal  modes 
of  a  dimer  and  the  translational  modes  in  dimer-monomer  collisions,  was  found  to  be  an  important  factor  that 
influences  the  magnitude  of  the  equilibrium  constant  Keq.  This  may  be  explained  as  follows.  The  dimers  are 
formed  after  three-body  collisions,  and  typically  have  internal  energies  smaller  than  the  evaporation  energy 
after  those  collisions.  In  argon,  the  evaporation  energy  for  a  dimer  is  relatively  small  compared  to  the  typical 
total  collision  energy  for  all  temperatures  under  consideration  (E^/k  «  140  K).  That  means  that  most  of 
the  dimers  will  have  their  internal  energy  in  excess  of  the  evaporation  energy  just  after  one  or  two  collisions 
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Figure  7.  Argon  dimer  formation  rate  as  a  function  of  gas  temperature. 


Figure  8.  Argon  dimer  formation  equilibrium  constant  as  a  function  of  gas  temperat 
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with  monomers.  The  lifetime  of  the  dimers  whose  internal  energy  is  larger  than  the  evaporation  energy 
is  very  short,  on  the  order  of  a  picosecond.  This  results  in  the  dimer-monomer  energy  transfer  being  the 
main  process  that  leads  to  quick  dimer  dissociation.  Note  that  the  value  of  Z  has  negligible  impact  on  the 
dimer  formation  rates,  and  only  the  evaporation  rates  are  affected.  As  a  result,  in  the  range  of  temperatures 
considered  in  this  work,  the  equilibrium  constant  for  argon  was  found  to  be  nearly  proportional  to  Z~l . 

The  Z  dependence  of  the  equilibrium  constant  is  weaker  for  water  molecule  condensation.  In  this  case, 
the  evaporation  energy  of  a  dimer  is  much  larger  than  the  translational  energy  of  colliding  molecules  and 
dimers  (the  reduced  evaporation  energy  E^/k  ~  3,500  K,  compared  to  gas  temperatures  on  the  order  of 
300  K).  The  high  value  of  the  evaporation  energy  results  in  longer  lifetimes  of  dimers,  since  much  more 
collisions  are  necessary  to  transfer  enough  energy  from  the  translational  modes  to  the  internal  modes  of  a 
dimer.  The  dependence  of  Keq  on  Z  is  therefore  much  weaker  for  water  than  for  argon.  In  a  250  K  thermal 
bath,  Keq  was  found  to  decreases  by  only  about  a  factor  of  two  when  Z  decreases  from  10  to  1. 

Comparison  of  the  equilibrium  constant  obtained  with  the  present  model  and  two  different  datasets,  with 
the  theoretical  results  of  Refs.  43,66,  where  a  flexible  potential  energy  surface  fitted  to  spectroscopical  data 
was  used,  is  shown  in  Fig.  9.  Note  that  there  are  a  number  of  theoretical  predictions  of  the  equilibrium  rate  of 
water  dimerization,  and  they  differ  by  at  least  a  factor  of  three  in  the  range  of  temperatures  considered  in  this 
work.  Reference  43  was  chosen  for  comparison  as  the  most  sophisticated  and  one  of  the  most  recent  ones.  The 
results  for  Dataset  1  show  that  the  calculated  equilirbium  constant  agrees  well  with  the  theorecal  prediction 
for  temperatures  160  K  and  250  K,  the  range  that  is  expected  to  be  very  important  in  terms  of  dimer 
formation  and  nucleation  in  plume  expansions.  The  calculated  values  for  higher  temperatures  significantly 
overpredict  the  theoretical  curve.  For  dataset  2,  the  results  are  close  to  Refs.  43,  66  for  temperatures 
between  200  K  and  350  K.  For  lower  temperatures,  the  calculated  values  are  noticeably  smaller.  The 
difference  between  Satasets  1  and  2  is  within  a  factor  of  two  for  the  entire  range  of  temperatures  under 
consideration,  which  indicates  that  the  change  in  the  binding  energy  and  heat  capacity,  that  are  assumed 
to  be  temperature  independent,  do  not  significantly  change  the  slope  of  the  equilibrium  constant  decreasing 
with  increasing  temperature.  This  general  trend  of  weaker  slope  may  be  related  to  a  number  of  factors, 
such  as  temperature  dependence  of  actual  heat  capacity  and  evaporation  energy,  as  well  as  the  approximate 
after-reaction  energy  redistribution  used  in  this  work  (recall  that  the  Larsen-Borgnakke  model  was  used  for 
the  energy  redistribution). 


Figure  9.  Water  equilibrium  constant  as  a  function  of  gas  temperature. 
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V.  Orifice  Expansion  Flow 


An  experimental  study67  of  dimer  formation  in  supersonic  water  vapor  molecular  beams  expanding  from 
a  circular  orifice  is  chosen  as  the  basis  for  comparison  and  analysis  of  the  present  algorithm  and  parameters 
of  the  numerical  model.  The  flow  of  water  vapor  through  a  123  /mum  diameter  orifice  is  modeled  for  a 
chamber  pressure  ranging  from  30  torr  to  300  torr,  and  two  reservoir  temperatures  of  373  K  and  493  K.  In  the 
DSMC  simulations,  the  subsonic  boundary  conditions  were  set  far  enough  from  the  orifice  plane  to  avoid  any 
impact  of  their  location.  Fully  diffuse  accommodation  on  the  wall  was  assumed.  To  avoid  slow  convergence 
associated  with  the  subsonic  part  of  the  domain,  the  calculations  included  to  steps.  First,  the  flow  was 
calculated  in  the  free  molecular  regime.  Then,  the  simulated  molecules  from  the  first  step  wer  utilized  as 
the  initial  condition  for  the  second,  high  pressure,  step.  The  typical  numbers  of  simulated  molecules  and 
collision  cells  were  20  million  and  2  million,  respectively.  Dataset  1  is  used  in  all  presented  results  except 
those  explicitely  indicated. 

The  gas  temperature  fields  and  the  cluster  number  density  fields  (all  cluster  sizes  included)  are  presented 
in  Fig.  10  for  the  smallest  and  highest  pressures  under  consideration.  While  the  gas  temperature  fields  are 
qualitatively  similar,  there  is  some  quantitative  difference  between  30  torr  and  300  torr  cases  that  is  expected 
to  be  substantial  in  terms  of  vapor  condensation.  In  the  near  field  of  the  orifice,  the  coreflow  temperatures 
are  significantly  lower  for  the  low  pressure  case.  It  is  closely  related  to  the  gas  mean  free  path  in  this  region. 
As  density  decreases  sharply  in  the  expansion,  there  is  relatively  few  collisions  in  the  plume  for  the  30  torr 
case.  This  results  in  larger  deflection  from  the  symmetry  axis,  and  as  a  result  lower  gas  temperatures  in  the 
coreflow.  Further  downstream,  larger  collision  frequency  results  in  lower  freezing  temperatures;  at  distances 
larger  than  1  mm  from  the  orifice  plane  the  gas  temperature  at  300  torr  is  lower  than  at  30  torr.  For  both 
pressures,  there  is  a  significant  number  of  clusters  in  the  plenum.  Generally,  their  number  corresponds  to 
the  equilibrium  constant  at  a  given  temperature  and  density;  it  is  two  orders  of  magnitude  larger  for  the 
higher  pressure  case.  The  descrease  in  the  plume  is  more  pronounced  for  this  case. 


Figure  10.  Impact  of  chamber  pressure  on  gas  translational  temperature  (K),  left,  and  cluster  number  density 
(m~3),  right.  Top  halves,  po  =  30  torr,  bottom  halves,  po  =  300  torr. 


This  decrease  is  illustrated  in  Fig.  11,  where  the  cluster  mole  fraction  normalized  by  its  chamber  value 
is  presented  along  with  the  water  vapor  temperature  for  two  plenum  pressures.  It  is  clearly  seen  that  the 
relative  increase  in  the  cluster  mole  fraction  is  more  significant  for  the  lower  pressure.  The  reason  for  this  is 
that  even  though  the  relative  gas  density  in  the  coreflow  is  higher  for  300  torr,  which  results  in  more  clusters 
produced,  it  is  not  compensated  by  significantly  lower  temperatures  in  the  expansion  region  for  this  case. 
For  example,  100  /mum  downstream  from  the  orifice  plane  the  temperature  for  300  torr  case  is  about  35K 
lower,  which  translates  to  over  five  times  lower  dimer  equilibrium  constant.  Most  of  the  condensation  occurs 
at  temperatures  between  300  K  and  150  K.  Note  that  at  about  200  /mum  from  the  orifice  plane  the  cluster 
mole  reaches  a  plateau,  at  which  point  any  further  increase  in  the  mole  fraction  along  the  symmetry  axis  is 
primarily  due  to  gasdynamic  reasons  (heavier  particles  tend  to  stay  in  the  coreflow),  and  not  as  much  due 
to  the  condensation-evaporation  mechanism.  For  a  higher  stagnation  temperature  of  493  K,  the  dimer  mole 
fraction  start  to  increase  at  about  450  K,  and  increases  until  the  gas  is  colled  down  to  about  150  K,  at  which 
point  further  increase  is  hampered  by  relatively  low  collision  frequencies. 

While  the  knowledge  of  total  cluster  mole  fractions  may  be  important  in  many  cases  and  applications, 
it  is  also  important  to  know  the  size  distribution  of  these  clusters.  Generally,  it  may  change  from  nearly 
Gaussian  at  low  pressures  to  bi-model  at  higher  pressures.33  The  terminal  size  distribution  of  water  clusters 
for  the  cases  under  consideration  is  shown  in  Fig.  12.  Wvwn  though  it  is  difficult  to  draw  any  conclusions 
for  the  two  smaller  pressures,  the  size  distribution  is  certainly  non-Gaussian  for  the  two  larger  pressures. 
The  upper  levels  are  clearly  underpopulated  as  compared  to  a  Boltzmann-like  equilibrium  distribution.  The 
reason  maybe  fast  cooling  and  simultaneous  increase  in  the  gas  mean  free  path  in  the  expanding  plume 
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Figure  11.  Gas  temperature  and  cluster  mole  fraction  profiles  along  the  orifice  axis  for  a  stagnation  temperature 
of  373  K  (left)  and  493  K  (right). 


flow,  where  there  are  not  enough  collisions  to  produce  significant  number  of  clusters  larger  than  4-mers. 
Nevertherless,  note  that  at  300  torr  the  relative  population  of  dimers  is  less  than  50%,  the  rest  is  larger 
clusters.  The  dimer  population  increases  to  over  60%  for  200  torr  and  to  about  85%  for  100  torr.  For  493  K 
cases  (not  shown  here),  the  relative  population  of  of  trimers  is  9%  for  300  torr,  2.5%  for  200  torr,  and  less 
than  1%  for  100  torr.  The  number  of  larger  clusters  is  negligible  for  this  temperature. 

Comparison  of  the  numerical  results  for  the  terminal  mole  fraction  of  dimers  near  the  symmetry  axis 
with  experimental  results67  is  presented  in  Fig.  13  (right)  for  two  considered  stagnation  temperatures.  For 
both  temperatures,  there  is  a  reasonable  agreement  for  lower  pressures  where  the  dimers  comprise  well  over 
90%  od  the  total  number  of  clusters.  For  higher  pressures,  where  the  important  processes  include  not  only 
the  dimer  formation  rate,  but  also  cluster  coalescence  and  monomer  sticking,  the  computational  results 
significantly  underpredict  the  data.  Additional  computations  conducted  for  373  K  with  the  coalescence 
probability  reduced  from  the  baseline  of  1  to  0.2,  show  the  terminal  dimer  mole  fraction  increase  to  0.02  at 
200  torr  and  0.021  at  300  torr.  Although  it  still  underpredicts  the  data,  it  is  over  a  factor  of  two  higher 
than  the  baseline  values.  The  calculations  were  also  performed  for  Dataset  1  (not  shown  here);  the  results 
are  generally  about  25%  lower  than  for  Dataset  2.  For  example,  for  100  torr  at  373  K  it  is  0.006  instead 
of  0.0083  for  Dataset  2.  This  was  expected  since  for  Dataset  1  the  dimer  formation  equilirbium  constant 
is  significantly  higher  at  low  temperatures.  There  may  be  a  number  of  possible  reasons  for  the  difference 
between  the  experimental  and  numerical  data  which  root  both  in  computational  and  measurement  procedures 
and  uncertainties.  Practically  all  parameters  used  in  the  present  model  and  listed  in  the  previous  section 
have  significant  error  bars.  Experimental  data  also  have  error  bars  related  to  the  assumed  cross  sections, 
intrusions  introduced  by  the  ionization  precedure.  The  actual  geometry  of  the  orifice,  such  as  the  orifice 
thickness,  as  well  as  orifice  surface  temperature  distribution,  are  also  unclear.  The  present  model  better 
captures  the  decrease  in  the  mole  fraction  due  to  the  stagnation  temperature  increase  than  the  model10  (see 
Fig.  13,  left),  but  again,  it  is  less  pronounced  than  in  the  experiment. 

VI.  Conclusions 

The  homogeneous  condensation  model  for  the  DSMC  method  is  considered  in  this  work  for  two  gases, 
argon  and  water.  Important  parameters  of  the  model,  such  as  the  binding  energy  and  the  heat  capacity  as 
function  of  cluster  size,  are  refined  as  compared  to  the  previous  work.23  For  argon,  the  present  model  was 
found  to  match  theoretical  dimer  formation  rate  and  equilirbium  constant  in  the  wide  range  of  temperautres. 
In  this  case,  a  temperature-dependent  inelastic  collision  number  for  monomer-cluster  interactions  was  intro¬ 
duced.  For  water,  two  datasets  were  used,  one  based  on  quantum  calculations  of  heat  capacity  and  binding 
energy  available  in  the  literature,  and  the  other  based  on  the  present  molecular  dynamics  computations. 
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Figure  12.  Terminal  cluster  size  distribution  function  for  different  pressures  at  a  stagnation  temperature  of 
373  K. 


For  both  datasets,  the  computed  equilibrium  constant  as  a  function  of  temperature  has  lower  slope  than 
the  available  theoretical  prediction.  Comparison  of  water  dimer  formation  in  a  plume  expanding  through 
a  circular  orifice,  with  experimental  data  shows  good  agreement  at  lower  pressures  and  significantly  lower 
numerical  values  for  higher  pressures.  Further  analysis  is  needed  to  clarify  the  reasons  for  such  a  difference. 
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